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ABSTRACT 

New paramcterizations for the spectral dissipation of wind-generated waves are proposed. The rates 
of dissipation have no predetermined spectral shapes and are functions of the wave spectrum, in a 
way consistent with observation of wave breaking and swell dissipation properties. Namely, swell 
dissipation is nonlinear and proportional to the swell steepness, and wave breaking only affects 
spectral components such that the non-dimensional spectrum exceeds the threshold at which waves 
are observed to start breaking. An additional source of short wave dissipation due to long wave 
breaking is introduced, together with a reduction of wind-wave generation term for short waves, 
otherwise taken from Janssen (J. Phys. Occanogr. 1991). These parameterizations are combined 
and calibrated with the Discrete Interaction Approximation of Hasselmann et al. (J. Phys. Oceangr. 
1985) for the nonlinear interactions. Parameters are adjusted to reproduce observed shapes of 
directional wave spectra, and the variability of spectral moments with wind speed and wave height. 
The wave energy balance is verified in a wide range of conditions and scales, from the global ocean 
to coastal settings. Wave height, peak and mean periods, and spectral data are validated using in 
situ and remote sensing data. Some systematic defects arc still present, but the parameterizations 
probably yield the most accurate overall estimate of wave parameters to date. Perspectives for 
further improvement are also given. 



1. Introduction 

a. On phase- averaged models 

Spectral wave modelling has been performed for the last 
50 years, using the wave energy balance equation (Gelci 
et al. 1957). This model for the evolution of spectral den- 



sities of the surface elevation variance F distributed over 
frequencies / and directions 9 can be put in the form 



dF(f,( 
dt 



= W/, 8)+Sm(f, 9)+S oc (f, 6)+S bt (f, 0), (1) 
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where the Lagrangian derivative is the rate of change of 
the spectral density when following a wave packet at its 
group speed in physical and spectral space. This spectral 
advection particularly includes changes in direction due to 
the Earth sphericity and refraction over varying topogra- 
phy (e.g. Munk and Traylor 1947; Magne et al. 2007) and 
currents, and changes in wavelength or period in similar 
conditions (Barber 1949). 

The source functions on the right hand side are sepa- 
rated into an atmospheric source function S atm , a nonlin- 
ear scattering term S n \, an ocean source S OCl and a bottom 
source Sbt- 

This separation, like any other, is largely arbitrary. For 
example, waves that break are highly nonlinear and thus 
the effect of breaking waves that is contained in S oc is in- 
trisically related to the non-linear evolution term contained 
in 5 n i. Yet, compared to the usual separation of deep-water 
evolution into wind input, non-linear interactions, and dis- 
sipation, it has the benefit of identifying where the energy 
and momentum is going to or coming from, which is a nec- 
essary feature when ocean waves are used to drive or are 
coupled with atmospheric or ocean circulation models (e.g. 
Janssen et al. 2004; Ardhuin et al. 2008b). 

<Satm, which gives the flux of energy from the atmo- 
spheric non-wave motion to the wave motion, is the sum 
of a wave generation term S- ln and a wind generation term 
S ut (often referred to as "negative wind input" , i.e. a 
wind output). The nonlinear scattering term 

^ni repre- 
sents all processes that lead to an exchange of wave energy 
and momentum between the different spectral components. 
In deep and intermediate water depth, this is dominated 
by cubic interactions between quadruplets of wave trains, 
while quadratic nonlinearities play an important role in 
shallow water (e.g. WISE Group 2007). The ocean source 
<S oc may accomodate wave-current interactions 1 and inter- 
actions of surface and internal waves, but it will be here 
restricted to wave breaking and wave-turbulence interac- 
tions. 

The basic principle underlying equation (1) is that waves 
essentially propagate as a superposition of almost linear 
wave groups that evolve on longer time scales as a result of 
weak-in-the-mean processes (e.g. Komen et al. 1994). Re- 
cent reviews have questioned the possibility of further im- 
proving numerical wave models without changing this basic 
principle (Cavaleri 2006). Although this may be true in the 
long term, we demonstrate here that it is possible to im- 
prove model results significantly by including more physical 
features in the source term parametcrizations. The main 
advance proposed in the present paper is the adjustment 
of a shape-free dissipation function based on today's em- 



1 In the presence of variable current, the source of energy for the 
wave field, i.e. the work of the radiation stresses, is generally hidden 
when the energy balance is written as an action balance (e.g. Komen 
et al. 1994). 



pirical knowledge on the breaking of random waves (Ban- 
ner et al. 2000; Babanin et al. 2001) and the dissipation 
of swells over long distances (Ardhuin et al. 2009b). The 
present formulations are still semi-empirical, in the sense 
that they are not based on a detailed physical model of dis- 
sipation processes, but they demonstrate that progress is 
possible. This effort opens the way for completely physical 
parametcrizations (e.g. Filipot et al. 2010) that will even- 
tually provide new applications for wave models, such as 
the estimation of statistical parameters for breaking waves, 
including whitecap coverage and foam thickness. Other ef- 
forts, less empirical in nature, are also under way to arrive 
at better parameterizations (e.g. Banner and Morison 2006; 
Babanin et al. 2007; Tsagareli 2008), but they yet have to 
produce a practical alternative for wave forecasting and 
hindcasting. 

b. Shortcomings of existing parameterizations 

All wave dissipation parameterizations up to the work 
of van der Westhuysen et al. (2007) had no quantitative re- 
lationship with observed features of wave dissipation, and 
the parameterizations were generally used as set of tuning 
knobs to close the wave energy balance. The parameteri- 
zation of the form proposed by Komen et al. (1984) have 
produced a family loosely justified by the so-called 'random 
pulse' theory of (Hasselmann 1974). These take a generic 
form 



Soc (f,6) = C ds g°- 5 k$- 5 Ht 



5l Y + 52 \k 



(2) 



in which Cd s is a negative constant, and k r is an energy- 
weighted mean wavenumber defined from the entire spec- 
trum, and H s is the significant wave height. In the early 
and latest parameterizations, the following definition was 
used 
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k r E(f,e)dfde 



(3) 



where r is a chosen real constant, typically r = —0.5 or 
r = 0.5. 

These parameterizations are still widely used in spite of 
inconsistencies in the underlying theory. Indeed, if white- 
caps do act as random pressure pulses, their average work 
on the underlying waves only occurs because of a phase cor- 
relation between the vertical orbital velocity field and the 
moving whitecap position, which travels with the breaking 
wave. In reality the horizontal shear is likely the dominant 
mechanism (Longuet-Higgins and Turner 1974), but the 
question of correlation remains the same. For any given 
whitecap, such a correlation cannot exist for all spectral 
wave components: a whitecap that travels with one wave 
leads to the dissipation of spectral wave components that 
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propagate in similar directions, with comparable phase ve- 
locities. However, whitecaps moving in one direction will 
give (on average) a zero correlation for waves propagating 
in the opposite direction because the position of the crests 
of these opposing waves are completely random with re- 
spect to the whitccap position. As a result, not all wave 
components arc dissipated by a given whitecap (others 
should even be generated), and the dissipation function 
cannot take the spectral form given by Komen et al. (1984). 

A strict interpretation of the pressure pulse model gives 
a zero dissipation for swells in the open ocean because the 
swell wave phases are uncorrelated to those of the shorter 
breaking waves. There is only a negligible dissipation due 
to short wave modulations by swells and preferential break- 
ing on the swell crests (Phillips 1963; Hasselmann 1971; 
Ardhuin and Jenkins 2005). Still, the Komen ct al. (1984) 
type dissipation terms are applied to the entire spectrum, 
including swells, without any physical justification. 

In spite of its successful use for the estimation of the 
significant wave height H s and peak period T p , these fixed- 
shape dissipation functions, from Komen et al. (1984) up 
to Bidlot et al. (2007a), have built-in defects. Most con- 
spicuous is the spurious amplification of wind sea growth 
in the presence of swell (e.g. van Vledder and Hurdle 2002), 
which is contrary to all observations (Dobson et al. 1989; 
Violante-Carvalho ct al. 2004; Ardhuin et al. 2007). As- 
sociated wih that defect also comes an underestimation of 
the energy level in the inertial range, making these wave 
models ill-suited for remote sensing applications, as will be 
exposed below. 

Also, these paramctcrizations typically give a decreas- 
ing dissipation of swell with increasing swell steepness, con- 
trary to all observations from Darbyshire (1958) to Ardhuin 
et al. (2009b). This effect is easily seen by taking a sea state 
composed of a swell and wind sea of energy E\ and E 2 and 
mean wavenumbers k\ and k 2 , respectively, with k 2 > k\. 
The overall mean wavenumber is 

k r = [(klEh + k r 2 E 2 ) /[Ex + E 2 )] 1/r . (4) 

Equation (2) gives a dissipation that is proportional to 
k^ 5 (Ei + E 2 ) in the low frequency limit. Now, if we keep 
k\ , &2 an E 2 constant and only increase the swell energy E\ , 
the relative change in dissipation is, according to (2), pro- 
portional to x = 3.5[(fci/fc r ) r — l]/r+2. For r = 0.5, as used 
by (Bidlot ct al. 2005, hereinafter BAJ), and kx/k r < 0.51, 
x is negative (i.e. the dissipation decreases with increasing 
swell energy). For equal energy in sea and swell, this occurs 
when k\jk 2 < 0.3, which is generally the case with sea and 
swell in the ocean. This erroneous decrease of swell dissi- 
pation with increasing swell steepness is reduced when the 
model frequency range is limited to maximum frequency of 
0.4 Hz, in which case the lowest winds (less than 5 m/s) are 
unable to produce a realistic wind sea level, hence limiting 
the value of k r to relatively small values. 



An alternative and widely used formulation has been 
proposed by Tolman and Chalikov (1996), and some of its 
features are worth noting. It combines two distinct dissipa- 
tion formulations for high and low frequencies, with a tran- 
sition at two times the wind sea peak frequency Whereas 
Jansscn ct al. (1994) introduced the use of two terms, k and 
k 2 in eq. (2), in order to match the very different balances 
in high and low frequency parts of the spectrum, they still 
had a common fixed coefficient, Cd s g°' 5 kf 5 H^. In Tolman 
and Chalikov (1996) these two dissipation terms are com- 
pletely distinct, the low frequency part being linear in the 
spectrum and proportional to wind friction velocity it*, the 
high frequency part is also linear and proportional to u 2 . In 
this formulation the frequency dependence of the two terms 
is also prescribed. Tolman and Chalikov (1996) further in- 
cluded swell attenuation by the wind, based on numerical 
simulations of the airflow above waves (Chalikov and Bclc- 
vich 1993), here noted S ou t- At relatively short fetches, 
these source terms are typically a factor 2 to 3 smaller than 
those of Janssen ct al. (1994), which was found to produce 
important biases in wave growth and wave directions at 
short fetch (Ardhuin et al. 2007). Another set of param- 
eterizations was proposed by Makin and Stam (2003). It 
is appropriate for high winds conditions but does not pro- 
duce accurate results in moderate sea states (Lefevre et al. 
2004). Finally, among the many formulations proposed we 
may cite one by Polnikov and Inocentini (2008), but its 
accuracy appears generally less than with the model pre- 
sented here, in particular for mean periods. 

Based on observations of large wave height gradients in 
rapidly varying currents, Phillips (1984) proposed a dissi- 
pation rate proportional to the non-dimensional spectrum 
B, also termed 'saturation spectrum'. Banner et al. (2000) 
indeed found a correlation between the breaking probabil- 
ity of dominant waves and the saturation, when the latter is 
integrated over a finite frequency bandwidth and all direc- 
tions., with breaking occuring when B exceeds a threshold 
B r . Alves and Banner (2003) proposed to define the dissi- 
pation S oc by B / B r to some power, multiplied by a Komcn- 
type dissipation term. Although this approach avoided the 
investigation of the dissipation of non-breaking waves, it 
imported all the above mentionned defects of that parame- 
terization. Further, these authors used a value for B r that 
is much higher than suggested by observations, which tends 
to disconnect the parameterization from the observed ef- 
fects (Babanin and van der Westhuysen 2008). 

The use of a saturation parameter was taken up again 
by van der Westhuysen ct al. (2007), hereinafter WZB, 
who, Like Alves and Banner (2003), integrated the satura- 
tion spectrum over directions, giving 

B(f)= k 3 F(f,6')C g /(2n)d6'. (5) 
Jo 
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From this, they defined the source function 

- B(/)lP /2 



:.WZB 



(f,6) = -C^gk 



(6) 



where C is a positive constant, B r is a constant saturation 
threshold and and p is a coefficient that varies both with 
the wind friction velocity u* and the degree of saturation 
B(f)/B r with, in particular, p w for £(/) < 0.8B r . For 
non breaking waves, when p«0, the dissipation is too large 
by at least one order of magnitude, making the parameteri- 
zation unfit for oceanic scale applications with wave heights 
in the Atlantic undcrprcdictcd by about 50% (Ardhuin and 
Le Boyer 2006). In van der Westhuysen (2007) this was 
addressed by reverting back to Komen ct al. (1984) dis- 
sipation for non-breaking waves, but no solution for the 
dissipation of these spectral components was proposed. In 
WZB, the increase of p with the inverse wave age u*/C 
was designed to increase S oc at high frequency, which was 
needed to obtain a balance with the S'atm term in equation 
(1). This indicates that, besides the value of the saturation 
B r , other factors may be important, such as the direction- 
ality of the waves (Banner et al. 2002). Other observations 
clearly show that the breaking rate of high frequency waves 
is much higher for a given value of B, probably due to cu- 
mulative effects by which the longer waves are modifying 
the dissipation of shorter waves. 

Banner et al. (1989) and Melville et al. (2002) have 
shown how breaking waves suppress the short waves on the 
surface, and we will show here that a simple estimation of 
the dominant breaking rates based on the observations by 
Banner et al. (2000) suggests that this effect is dominant 
for wave frequencies above three times the windsea peak 
frequency. Young and Babanin (2006) arrived at the same 
conclusion from the examination of wave spectra, and pro- 
posed a parameterization for S oc that included a new term, 
the cumulative term, to represent theis effect. Yet, their 
estimate was derived for very strong wind-forcing condi- 
tions only. Further, their interpretation of the differences 
in parts of a wave record with breaking and non-breaking 
waves implies an underestimation of the dissipation rates 
because the breaking waves have already lost some energy 
when they are observed and the non-breaking waves are not 
going to break right after they have been observed. Also, 
since the spectra arc different, nonlinear interactions must 
be different, even on this relatively small time scale (e.g. 
Young and van Vledder 1993, figure 5), and the differences 
in spectra may not be the result of dissipation alone. 

Finally, the recent measurement of swell dissipation by 
Ardhuin et al. (2009a) has revealed that the dissipation of 
non-breaking waves is essentially a function of the wave 
steepness, and a very important process for ocean basins 
larger than 1000 km. Because of the differences between 
coastal and larger scale sea states (e.g. Long and Resio 
2007), it is paramount to verify the source function pa- 



rameterizations at all scales, in order to provide a robust 
and comprehensive parameterization of wave dissipation. 

c. A new set of parameterizations and adjustments to get ade- 
quate balances 

It is thus time to combine the existing knowledge on the 
dissipation of breaking and non-breaking waves to provide 
an improved parameterization for the dissipation of waves. 
Our objective is to provide a robust parameterization that 
improves existing wave models. For this we will use the 
parameterization by BAJ as a benchmark, because it was 
shown to provide the best forecasts on global scales (Bidlot 
et al. 2007b) before the advent of the parameterizations 
presented here. BAJ is also fairly close to the widely used 
"WAM-Cycle 4" parameterization by Janssen and others 
(Komen et al. 1994). 

We will first present a general form of the dissipation 
terms based on observed wave dissipation features. The 
degrees of freedom in the parameterization are then used 
to adjust the model result. In particular we adujst a cu- 
mulative breaking effect and a wind sheltering effect that, 
respectively, dissipates and reduces the wind input to short 
waves as a function of longer waves characteristics. A com- 
prehensive validation of wave parameters is then presented 
using field experiments and a one year hindcast of waves 
at the global and regional scale, in which all possible wave 
measurements are considered, with significant wave heights 
ranging from to 17 m. The model is further validated 
with independent data at regional and global scales. 

Tests and verification in the presence of currents, and 
using a more realistic parameterizations of wave-wave in- 
teractions will be presented in parts II and III. These may 
also include some replacement of the arbitrary choices made 
here in the details of the dissipation parameters, with physically- 
motivated expressions. 

2. Parameterizations 

Several results will be presented, obtained by a nu- 
merical integration of the energy balance. Because nu- 
merical choices can have important effects (e.g. Tolman 
1992; Hargreaves and Annan 2000), a few details should 
be given. All calculations are performed with the WAVE- 
WATCH III™ modelling framework (Tolman 2008, 2009), 
hereinafter WWATCH, using the third order spatial and 
spectral advection scheme, and including modifications of 
the source terms described here. In all cases ran with 
WWATCH, the source terms are integrated with the fully 
implicit scheme of Hargreaves and Annan (2000), combined 
with the adaptative time step and limiter method of Tol- 
man (2002a), in which a minimum time step of 10 s is 
used, so that the limiter on wind-wave growth is almost 
never activated. The diagnostic tail, proportional to / 
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is only imposed at a cut-off frequency f c set to 

fc = /fm/tti- (7) 

Here we take /fm = 10 and define the mean frequency as 
/m=l/7m0,i- Hence f c is generally above the maximum 
model frequency that we fixed at 0.72 Hz, and the high 
frequency tail is left to evolve freely Some comparison tests 
are also done with other parameterizations using a lower 
value of f c , typically set at 2.5 f m (Bidlot et al. 2007a). 
In such calculations, although the net source term may 
be non-zero at frequencies above f c , there is no spectral 
evolution due to the imposed tail. 

a. Nonlinear wave wave interactions 

All the results discussed and presented in this section 
are obtained with the Discrete Interaction Approximation 
of Hassclmann et al. (1985). The coupling coefficient that 
gives the magnitude of the interactions is G n \. Based on 
comparisons with exact calculations, Komen et al. (1984) 
adjusted the value of C n i to 2.78 x 10 7 , which is the value 
used by Bidlot et al. (2005). Here this constant will be al- 
lowed to vary slightly. This parameterization is well known 
for its shortcomings (Banner and Young 1994), and the 
adjustment of other parameters probably compensates for 
some of these errors. This matter will be fully discussed in 
Part III. 

b. Swell dissipation 

Observations of swell dissipation are consistent with the 
effect of friction at the air-sea interface (Ardhuin et al. 
2009a), resulting in a flux of momentum from the wave 
field to the wind (Harris 1966). We thus write the swell 
dissipation as a negative contribution S ou t which is added 
to Si n to make the wind- wave source term S^tm- 

Using the method of Collard et al. (2009), a systematic 
analysis of swell observations by Ardhuin et al. (2009a) 
showed that the swell dissipation is non-linear, possibly 
related to a laminar-to-turbulent transition of the oscilla- 
tory boundary layer over swells. Defining the boundary 
Reynolds number Re= 4u or ba rb/f / a, where w or b and a or b 
are the significant surface orbital velocity and displacement 
amplitudes, and v a is the air viscosity, we take, for Re less 
than a critical value Re c 

S ont (f,8) = -C dsv ^{2kV2^}F(f,9), (8) 

Pw L J 

where the constant Cdsv is equal to 1 in Dore (1978) 's lam- 
inar theory. 

When the boundary layer is expected to be turbulent, 
for Rc>Rc c , we take 

Sout {f, 9) = - ^ { 16 f e a 2 u olh /g} F (/, 9) . (9) 

Pw 



A few tests have indicated that a threshold Re c = 2 x 
10 5 m/H s provides reasonable result, although it may be 
a also be a function of the wind speed, and we have no 
explanation for the dependence on H s . A constant thresh- 
old close to 2 x 10 5 provides similar results. Here we shall 
use Cdsv = 1-2, but the results are not too sensitive to the 
exact value. 

The parameterization of the turbulent boundary layer 
is more problematic. Without direct measurements in the 
boundary layer, there is ample room for speculations. From 
the analogy with an oscillatory boundary layer over a fixed 
bottom (Jensen et al. 1989), the values of f e inferred from 
the swell observations, in the range 0.004 to 0.013 (Ardhuin 
et al. 2009b), correspond to a surface with a very small 
roughness. Because, we also expect the wind to influence 
/ e , the parameterization form includes adjustable effects 
of wind speed on the roughness, and an explicit correction 
of f e . This latter correction takes the form of a Taylor 
expansion to first order in u^/u^. 

fe = *1 \fe,GM + [M + *2 COs(0 - 6 U )] —\ , (10) 
I "orb J 

where f e ,GM is the friction factor given by Grant and Mad- 
sen's (1979) theory for rough oscillatory boundary layers 
without a mean flow. Adequate swell dissipation is ob- 
tained with constant values of f e in the range 0.004 to 
0.007, but these do not necessarily produce the best re- 
sults when comparing wave heights to observations. Based 
on the simple idea that most of the air-sea momentum flux 
is supported by the pressure-slope correlations that give 
rise to the wave field (Donelan 1998; Peirson and Banner 
2003), we have set the surface roughness to 

4 = r z0 z (11) 

where r z o is here set to 0.04, of the roughness for the wind. 
We thus give the more generic equation (10) for f e , with 
fe,GM of the order of 0.003 for values of a or b / z' of the order 
of 2xl0 5 . 

The coefficients S2 and S3 of the O(it+/u or b) correc- 
tion have been adjusted to -0.018 and 0.015, respectively, 
the former negative value giving a stronger dissipation for 
swells opposed to winds, when cos (6* — 9 U ) < 0. This gives 
a range of values of f e consistent with the observations, 
and reasonable hindcasts of swell decay (Fig. 1), with a 
small underestimation of dissipation for steep swells. An 
increase of si from 0.8 to 1.1 produces negative biases on 
H s of the order of 30% at all oceanic buoys (40% for a par- 
tial wave height estimated from a spectrum restricted to 
periods around 15 s), so that the magnitude of the swell dis- 
sipation cannot be much larger than chosen here. Further 
discussion and validation of the swell dissipation is pro- 
vided by the global scale hindcasts in Section 4. 
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Fig. 1. Comparison of modelled swell significant heights, 
following the propagation of the two swells shown by Ard- 
huin et al. (2009) with peak periods of 15 s and high and 
low dissipation rates. Known biases in the level 2 data have 
been corrected following Collard et al. (2009). 



c. Wave breaking 

Observations show that waves break when the orbital 
velocity at their crest U c comes close to the phase speed 
C, with a ratio U c /C > 0.8 for random waves (Tulin and 
Landrini 2001; Stansell and MacFarlane 2002; Wu and Nepf 
2002). It is nevertheless difficult to parameterize the break- 
ing of random waves, since the only available quantity here 
is the spectral density. This density can be related to the 
orbital velocity variance in a narrow frequency band. This 
question is addressed in detail by Filipot et al. (2010). Yet, 
a proper threshold has to be defined for this quantity, and 
the spectral rate of energy loss associated to breaking has 
to be defined. Also, breaking is intricately related to the 
complex non-linear evolution of the waves (e.g. Banner and 
Peirson 2007). 

These difficulties will be ignored here. We shall param- 
eterize the spectral dissipation rate directly from the wave 
spectrum, in a way similar to WZB. Essentially we distigu- 
ish between spontaneous and induced breaking, the lat- 
ter being caused by large scale breakers overtaking shorter 
waves, and causing them to be dissipated. For the spon- 
taneous breaking we parameterize the dissipation rate di- 
rectly from the spectrum, without the intermediate step of 



estimating a breaking probability. 

We started from the simplest possible dissipation term 
formulated in terms of the direction-integrated spectral sat- 
uration B (/) given by eq. (5), with a realistic threshold 
B 0r = 1.2 x 10~ 3 corresponding to the onset of wave break- 
ing (Babanin and Young 2005) . This saturation parameter 
corresponds exactly to the a parameter defined by Phillips 
(1958). The value B = 8 x 10~ 3 , given by Phillips, cor- 
responds to a self-similar sea state in which waves of all 
scales have the same shape, limited by the breaking limit. 

This view of the sea state, however, ignores completely 
wave directionality. Early tests of parameterizations based 
on this definition of B indicated that the spectra were too 
narrow (Ardhuin and Le Boyer 2006). This effect could 
be due to many errors. Because Banner et al. (2002) in- 
troduced a directional width in their saturation to explain 
some of the variability in observed breaking probabilities, 
we similarly modify the definition of B. Expecting also 
to have different dissipation rates in different directions, 
we define a saturation B' that would correspond, in deep 
water, to a normalized velocity variance projected in one 
direction (in the case Sb = 2), with a further restriction of 
the integration of directions controlled by Ag, 



B'(f,9) 



e+A e 



k 3 cos SB 



9-A e 



9')F(f,9')^d6', 

Z7T 



(12) 



Here we shall always use A s = 80°. As a result, a sea state 
with two systems of same energy but opposite direction will 
typically produce much less dissipation than a sea state 
with all the energy radiated in the same direction. 

We finally define our dissipation term as the sum of 
the saturation-based term of Ardhuin et al. (2008a) and a 
cumulative breaking term 5bk ;CU , 



S oc (f,6) 



- ( a j 



5 d max{B(f)-B r ,0Y 

(1 - S d ) max{5' (/, 9) - B r , 0} 2 ] F(f, 6) 

Sbk,cu(/,0)+S turb (/,0). (13) 



where 



B(f) =max{B'(/,f3),(9e [0,2tt[} 



(14) 



The combination of an isotropic part (the term that 
multiplies S^) and a direction-dependent part (the term 
with 1 — Sd) was intended to allow some control of the 
directional spread in resulting spectra. This aspect is illus- 
trated in figure 2 with a hindcast of the November 3 1999 
case during the Shoaling Waves Experiment (Ardhuin et al. 
2007). Clearly, the isotropic saturation in the TEST442 
dissipation (with the original threshold B r = 0.0012) pro- 
duces very narrow spectra, even though it is known that 
the DIA parameterization for nonlinear interactions, used 
here, tends to broaden the spectra. The same behaviour 
is obtained with the isotropic parameterization by van dcr 
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Fig. 2. Wave spectra on 3 November 1999 at buoy X3 
(fetch 39 km, wind speed U w = 9.4 m s -1 ), averaged over 
the time window 1200-1700 EST, from observations and 
model runs, with different model parameterizations (sym- 
bols): BAJ stands for Bidlot et al. (2005). (a) Energy, (b) 
mean direction (c) directional spread. This figure is ana- 
logue to the figures 10 and 11 in Ardhuin et al. (2007), the 
model forcing and setting are identical. It was further ver- 
ified that halving the resolution from 1 km to 500 m does 
not affect the results. All parameters for BAJ, TEST441 
and TEST443 are listed in tables 3 and 4. Input parame- 
ters for TEST443 are identical to those for TEST441, and 
TEST442 differs from TEST441 only in its isotropic di- 
rect breaking term, given by sb = 0, Ag = 180°, and 
B r = 0.0012. It should be noted that the overall dissi- 
pation term in TEST443 is made anisotropic due to the 
cumulative effect, but this does not alter much the under- 
estimation of directional spread. 



Westhuysen et al. (2007) , as demonstrated by Ardhuin and 
Lc Boyer (2006). Further, using an isotropic dissipation at 
all frequencies yields an energy spectrum that decays faster 
towards high frequencies than the observed spectrum (Fig. 
2. a). On the contrary, a fully directional dissipation term 
(TEST443 with 6 a = 0) gives a better fit for all parameters. 
With sb = 2, we reduce B r to 0.0009, a threshold for the 
onset of breaking that is consistent with the observations 
of Banner et al. (2000) and Banner et al. (2002). 

The dissipation constant C^' was adjusted to 2.2 x 
10 -4 in order to give acceptable time-limited wave growth 
and reasonable directions in fetch-limited growth (Ard- 
huin et al. 2007). As noted in this previous work, similar 
growths of wave energy with fetch are possible with almost 
any magnitude of the wind input, but a reasonable mean 
direction in slanting fetch conditions selects the range of 
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Fig. 3. Fetch-limited growth of the windsea energy as 
a function of fetch on 3 November 1999, averaged over 
the time window 1200-1700 EST, from observations and 
model runs, with different model parameterizations (sym- 
bols): BAJ stands for Bidlot et al. (2005). This figure is 
analogue to the figure 8 in Ardhuin et al. (2007) , the model 
forcing and setting are identical. All parameters for BAJ, 
TEST441 and TEST443 arc listed in tables Al and A2. 
Input parameters for TEST443 are identical to those for 
TEST441, and TEST442 differs from TEST441 only in its 
isotropic direct breaking term, given by sb = 0, Ag = 180°, 
and B r = 0.0012. 



possible levels of input. Here the mean directions at the 
observed peak frequency are still biassed by about 25° to- 
wards the alongshore direction with the parameterizations 
proposed here (Fig. 2.b), which is still less that the 50° ob- 
tained with the weaker Tolman and Chalikov (1996) source 
terms (Ardhuin et al. 2007, figure 11). A relatively better 
fit is obtained with the BAJ parameterization. This is 
likely due to either the stronger wind input or the weaker 
dissipation at the peak. It is likely that both features of 
the BAJ parameterization are more realistic than what we 
propose here. 

Figure 3 shows the fetch-limited growth in wave energy 
of various parameterizations. We repeat here the sensi- 
tivity test to the presence of swell, already displayed in 
Ardhuin et al. (2007). Whereas the 1 m swell causes an 
unrealistic doubling of the wind sea energy at short fetch 
in the BAJ parameterization, the new parameterizations, 
just like the one by van der Westhuysen et al. (2007) are, 
by design, insensitive to swell (not shown). 

The dissipation 5turb due to wave-turbulence interac- 
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tions (Ardhuin and Jenkins 2006) is expected to be much 
weaker than all other terms and will be neglected here. 

Finally, following the analysis by Filipot et al. (2010), 
the threshold B r is corrected for shallow water, so that 
B' /B r in different water depths corresponds to the same 
ratio of the root mean square orbital velocity and phase 
speed. For periodic and irrotational waves, the orbital ve- 
locity increases much more rapidly than the wave height as 
it approaches the breaking limit. Further, due to nonlinear 
distortions in the wave profile in shallow water, the height 
can be twice as large as the height of linear waves with the 
same energy. In order to express a relevant threshold from 
the elevation variance, we consider the slope kH\ m (kD) of 
an hypothetical linear wave that has the same energy as 
the wave of maximum height. In deep water 2 , kH\ m (oo) m 
0.77, and for other water depths we thus correct B r by a 
factor (kH\\ n (kD) / Hy m (oo)) 2 . Using streamfunction the- 
ory (Dalrymplc 1974), a polynomial fit as a function of 
Y = tanh(fe-D) gives 

B' r = B r Y [M 4 Y 3 + M 3 Y 2 + M 2 Y + Mi] . (15) 

such that B' r = B r in deep water. The fitted constants 
are M 4 = 1.3286, M 3 = -2.5709 , M 2 = 1.9995 and 
Mi = 0.2428. Although this behaviour is consistent with 
the variation of the depth-limited breaking parameter 7 
derived empirically by Ruessink et al. (2003), the result- 
ing dissipation rate is not yet expected to produce realistic 
results for surf zones because no effort was made to verify 
this aspect. This is the topic of ongoing work, outside of 
the scope of the present paper. 

The cumulative breaking term Sbk,cu represents the smooth- 
ing of the surface by big breakers with celerity C that wipe 
out smaller waves of phase speed C. Due to uncertainties 
in the estimation of this effect in the observations of Young 
and Babanin (2006) , we use the theoretical model of Ard- 
huin et al. (2009b). Briefly, the relative velocity of the 
crests is the norm of the vector difference, Ac = |C — C'|, 
and the dissipation rate of short wave is simply the rate of 
passage of the large breaker over short waves, i.e. the inte- 
gral of AcA(C)dC, where A(C)dC is the length of break- 
ing crests per unit surface that have velocity components 
between C x and C x + dC x , and between C y and C y + dC y 
(Phillips 1985). Because there is no consensus on the form 
of A (Gemmrich et al. 2008), we prefer to link A to break- 
ing probabilities. Based on Banner et al. (2000, figure 
6, br = 22 (e — 0.055) ), and taking their saturation pa- 
rameter e to be of the order of 1.6-y/ £?'(/, 6), the breaking- 
probability of dominant waves waves is approximately 

P = 56.8 (max{y/ B'(f, 9) - ^W^Q}^ . (16) 

2 This value of the maximum equivalent linear height Hy ln = 
2\J1E, with E the elevation variance, is smaller than the usual value 
kH = 0.88 due to the correction for the nonlinear wave profile for 
which H > \/2E. 



However, because they used a zero-crossing analysis, for a 
given wave scale, there are many times when waves are not 
counted because the record is dominated by another scale: 
in their analysis there is only one wave at any given time. 
This tends to overestimate the breaking probability by a 
factor of 1.5 to 2 (Manasseh et al. 2006), compared to the 
present approach in which we consider that several waves 
(of different scales) may be present at the same place and 
time. We shall thus correct for this effect, simply dividing 
Phy 2. 

With this approach we define the spectral density of 
crest length (breaking or not) per unit surface Z(k) such 
that J l(k)dk x dk y is the total length of all crests per unit 
surface, with a crest being defined as a local maximum of 
the elevation in one horizontal direction. In the wavenum- 
ber vector spectral space we take 

Z(k) = l/(27r 2 fc) (17) 

which is equivalent to a constant in wavenumber-direction 
space l(k,9) = l/(27r 2 ). This number was obtained by 
considering an ocean surface full of unidirectional waves, 
with one crest for each wavelength 2ir/k for each spectral 
interval Ak = k, e.g. one crest corresponding to spectral 
components in the range 0.5 k to 1.5 k. Wc further double 
the potential number of crests to account for the direction- 
ality of the sea state. These two assumptions have not 
been verified and thus the resulting value of Z(k) is merely 
an adjustable order of magnitude. 

Thus the spectral density of breaking crest length per 
unit surface is A(k) = Z(k)P(k). Assuming that any break- 
ing wave instantly dissipates all the energy of all waves 
with frequencies higher by a factor r cu or more, then the 
cumulative dissipation rate is simply given by the rate at 
which these shorter waves are taken over by larger breaking 
waves, times the spectral density, namely 

Sbk,cu(/,0) = C m F(f,0) f A c A(k')dk', (18) 

Jf'<Tcnf 

where r cu defines the maximum ratio of the frequencies of 
long waves that will wipe out short waves. 

We now obtain A by extrapolating eq. (16) to higher 
frequencies, 

Sbk.cu (/, 9) = C CU F (/, 9) J rcu/ J 27r 

xmax ^B{f',9') - /B;,o| 2 ^dO'df, 

(19) 

We shall take r cu = 0.5, and C cu is a tuning coefficient 
expected to be a negative number of order 1, which also 
corrects for errors in the estimation of I. 

This generic form of the source terms produces markedly 
different balances for both mature and fully developped 
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Fig. 4. Academic test case over a uniform ocean with a 
uniform 10 m s -1 wind starting from rest, after 8 hours 
of integration, when C p /Uiq f=a 1. Source term balances 
given by the parameterization BAJ, and the paramctcr- 
izations proposed here with the successive introduction 
of the cumulative breaking and the wind sheltering ef- 
fects with the parameters C cu and s u . For BAJ, a di- 
agnostic f~ 5 tail is applied above 2.5 the mean frequency. 
In order to make the high frequency balance visible, the 
source terms are multiplied by the normalization function 
M(f) = p w C/(p a E(f)aUio). The result with the parame- 
terization of Tolman and Chalikov (1996) is also given for 
reference. 

seas. For mature seas, without cumulative effect, figure 
(4) shows that a balance is possible that gives roughly the 
same energy level and wind input term as the BAJ param- 
eterization, up to 0.4 Hz. However, the balance for higher 
frquencies produces energy levels decreasing slower than 
/~ 4 as the dissipation is too weak compared to the input, 
and thus the nonlinear energy flux is reversed, pumping 
energy from the tail to lower frequencies. 

The introduction of a strong cumulative term (TEST437) 
allows a balance at roughly the same energy level. How- 
ever, with the present formulation this will lead to a dissi- 
pation too strong at high frequency for higher winds. The 
introduction of the sheltering effect via the parameter s u 
(details in section d) is designed to get a balance with a 
weaker cumulative effect. 



The most important qualitative feature is the lack of 
a regular predefined shape for the normalized dissipation 
term S oc (f)/E(f). Whereas the shape given by 5k + (1 — 
S)k 2 is clearly visible in BAJ (with extremely high dissipa- 
tion rates if one considers high frequencies), and the low to 
high frequency dissipation transition at 2 f p is evident in 
TC, the shape of the new dissipation rates are completely 
dictated by the local spectral saturation level. This leads 
to a relatively narrow peak of dissipation right above the 
spectral peak, where saturation is strongest. 

This feature helps to produce the realistic spectral shapes 
near the peak, with a steeper low frequency side and a more 
gentle slope on the high frequency side, contrary to the 
backward facing spectra produced by BAJ and TC. How- 
ever, this localized strong relative dissipation, S oc (f)/E(f), 
is hard to reconcile with time and spatial scales of break- 
ing events, and thus probably exaggerated. Indeed, there 
should be no significant difference in relative dissipation 
among the spectral components that contribute to a break- 
ing wave crest, provided that they do not disperse signifi- 
cantly over the breaker life time, which is less than a wave 
period. Linear dispersion of waves with frequencies that 
differ by only 10% should only produce small relative phase 
shift. Thus, there is no physical reason why a breaking 
event would take much more energy, relatively speaking, 
from the spectral band 1.1 to 1.2 f p than from 1.2 to 1.3/ p . 
The factor 2 difference produced here in the relative dissi- 
pation rates is unrealistic. This strong relative dissipation 
at the peak (50% higher than in BAJ) is one important 
factor that gives a slower growth of the wave spectrum 
in TEST441 compared to BAJ. It is possible that the lo- 
calization of S oc at the peak compensates for the broader 
spectrum produced by the DIA compared to results with 
an exact non-linear interaction calculation. 

We now consider "fully developped" conditions, illus- 
trated by figure 5, corresponding to long durations with 
steady wind and infinite fetch. At low frequency, the non- 
linear swell damping term S ou t (the negative part of Si„) 
cancels about 30 to 50% of the nonlinear energy flux, so 
that the sea state grows only very slowly. As a result 
"full development" does not exist (the wave height keeps 
growing), but the resulting energy is still compatible with 
the observations of mature sea states (Alves et al. 2003). 
In contrast, the linear swell damping adjusted by Tolman 
(2002b) to produce reasonable swell heights in the trop- 
ics is much smaller than the non-linear energy flux to low 
frequencies, even with the reduced interaction coefficient 
proposed by Tolman and Chalikov (1996). A non-linear 
swell dissipation appears necessary to obtain both a realis- 
tic damping of observed swells and a satisfatory agreement 
with mature wind waves. Nonlincarity also bring within 
the same order of magnitude the decay scales estimated 
for short (Hogstrom et al. 2009) and very long swells (Ard- 
huin et al. 2009b). 
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Fig. 5. Same as figure (4) but after 48 h of integration 
and without normalization of the source terms. Source 
term balances given by the parameterization BAJ, and the 
parametcrizations proposed here with the successive intro- 
duction of the cumulative breaking and the wind sheltering 
effects with the parameters C cu and s u . For BAJ, a diag- 
nostic f~ 5 tail is applied above 2.5 the mean frequency. 



Both parametcrizations are physically very different from 
the parameterizations of the Komen et al. (1984) family, in- 
cluding Bidlot ct al. (2005). In these, the swell energy is 
lost to the ocean via whitccapping. Here we propose that 
this energy is lost to the atmosphere, with an associated 
momentum flux that drives the wave-driven wind observed 
in laboratories (Harris 1966) and for very weak winds at 
sea (Smedman et al. 2009). 

In the inertial range, a reasonable balance of all the 
source terms is obtained for C cu = —0.4 (figure 5). In this 
case, the spectrum approaches an f~ A shape, for which S n i 
goes to zero, whereas in BAJ it decreases even faster than 
/ _5 which makes S n i positive above 0.4 Hz. The behaviour 
of the high frequency tail is best seen when displayed in 
non-dimensional form, as done in figure 6. This shows the 
unrealistic high level of the tail without cumulative effect 
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Fig. 6. Values of the spectral saturation Bq for the cases 
presented in figure 5. 

nor modification of the wind input, and the equally unre- 
alistic low tail with the BAJ parameterization, especially 
when the tail is left to evolve freely. 

Adding the cumulative effect can be used to control 
the tail level, but this degree of freedom is not enough. 
Indeed, in strongly forced conditions the dominant waves 
break frequently, and a high cumulative effect, C cu = — 1, 
reduces the energy level in the tail below observed lev- 
els. This effect can be seen by considering satellite-derived 
mean square slopes (Fig. 8), or high moments of the fre- 
quency spectrum derived from buoy data (not shown but 
similar) . 

That effect can be mitigated by decreasing C cu or in- 
creasing r cu , so that dominant breaking waves will only 
wipe out much smaller waves. Instead, and because the 
wind to wave momentum flux was apparently too high in 
high winds, we chose to introduce one more degree of free- 
dom, allowing a reduction of the wind input at high fre- 
quency. 

d. Wind input 

The wind input parameterization is thus adapted from 
Jansscn (1991) and the following adjustments performed 
by Bidlot et al. (2005, 2007a). The full wind input source 
term reads 



Sin(f,e) 



silf(f,e) + 



Pa Pmax z y4 
Pw K 



x max{cos(0 - 6 U ),0} P aF (/, 0) , (20) 



where /3 max is a non-dimensional growth parameter (con- 
stant), k is von Karman's constant. In the present im- 
plementation the air/water density ratio is constant. The 
power of the cosine is taken constant with p = 2. We de- 
fine the effective wave age Z = log(/i) where p, is given 
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by Janssen (1991), and corrected for intermediate water 
depths, so that 

Z = log(fczi) + k/ [cos (6 - 9 U ) (u+/C + Za)] , (21) 

where z\ is a roughness length modified by the wave-supported 
stress t w , and z a is a wave age tuning parameter. z\ is im- 
plicitly defined by 



C/io = 

zo = 
z\ = 



u* . I z u 
— log — 

K \Zi 



mm <^ a -,z 0max 
9 



x/T- t w /t' 



(22) 

(23) 
(24) 



where z u is the height at which the wind speed is spec- 
ified, usually 10 meters. The maximum value of Zq was 
added to reduce the unrealistic wind stresses at high winds 
that are otherwise given by the standard parameteriza- 
tion. For example, zo,max = 0.0015 is equivalent to setting 
a maximum wind drag coefficient of 2.5 x 10 -3 . For the 
TEST441 parameterization, we have adjusted z a = 0.006 
and A nax = 1.52 (Fig. 7). 

An important part of the parameterization is the cal- 
culation of the wave-supported stress t w , which includes 
the resolved part of the spectrum, as well as the growth 
of an assumed f~ 5 diagnostic tail beyond the highest fre- 
quency. This parameterization is highly sensitive to the 
high frequency part of the spectrum since a high energy 
level there will lead to a larger value of z\ and it*, which 
gives a positive feedback and reinforces the energy levels. 

In order to allow a balance with the saturation-based 
dissipation, the wind input at high frequency is reduced by 
modifying the friction velocity u±. This correction also 
reduces the drag coefficient at high winds. Essentially, 
the wind input is reduced for high frequencies and high 
winds, loosely following Chen and Belcher (2000). This is 
performed by replacing in eq. (20) with a frequency- 
dependent u'^{f) defined by 



MY 



ul (cos 9 u ,sm6 u ) 

s in (f',e) 



o Jo 



2 71 



c 



(cos 9, sin 9) dfd0, 



(25) 



where the sheltering coefficient \s u \ ~ 1 can be used to tune 
the wind stresses for high winds, which would be largely 
overestimated for s u = 0. For s u > this sheltering is 
also applied within the diagnostic tail, which requires the 
estimation of a 3-dimensional look-up table for the high 
frequency stress. The shape of the new wind input is il- 
lustrated in figure 7 for fully developped seas. Clearly, 
for relatively young waves the energy levels at the spectral 
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Fig. 7. Incremental adjustements to the wind- wave in- 
teraction source term <S a tmj going from the BAJ form to 
the one used in TEST441. From one curve to the next, 
only one parameter is changed. The computations are per- 
formed for the same spectrum obtained by running the 
model from a calm sea for 8 hours with the BAJ param- 
eterization and a wind speed of 10 m s _1 . The reduction 
of z a from 0.011 to 0.006 strongly reduces the input for 
frequencies in the range 0.15 to 0.2 Hz, which is probably 
overestimated in BAJ when average (20%) levels of gusti- 
ness are considered: now the wind input goes to zero for 
/ = 0.13, which corresponds to C/Uio = 0.83, whereas 
it is still significant at that wave age in the BAJ parame- 
terization. As a result the much lower input level need a 
readjustment, performed here by increasing /3 max to 1.52. 
Yet, this high value of /3 max produces very high wind stress 
values and thus a very strong high frequency input. Adding 
the sheltering term s u = 1 allows a decent balance at high 
frequency. Finally, the addition of the air-sea friction term 
that gives swell dissipation produces a significant reduction 
of the input to the wind sea at / = 0.25 Hz. It is question- 
able whether this mechanism also applies in the presence 
of the critical layer for those waves. This matter clearly 
requires more theoretical and experimental investigation. 

peak are lower with s u = 0.4 (TEST441) than in other 
runs, this is largely due to a reduced feedback of the wave 
age on the wind stress via the t w /t term in eq. (24). 

3. Consequences of the source term shape 

We have already illustrated the effects of various param- 
eters on spectral shapes in academic time-limited and more 
realistic fetch-limited conditions. We now look at real sea 
states observed in the world ocean. Although wave spectra 
are difficult to compare to the few available observations, 
we have investigated the systematic variation of spectral 
moments 



fo 



f n E(f)df. 



(26) 
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with 7i = 2,3 and 4, and cut-off frequencies in the range 
0.2 to 0.4 Hz. Such moments are relevant to a variety of 
applications. Ardhuin et al. (2009b) investigated the third 
moment, which is proportional to the surface Stokes drift 
in deep water, and found that buoy data are very well 
represented by a simple function, which typically explains 
95% of the variance, 



m 3 (fc) 



5.9gU_ 



10 



x 10~ 



(2tt) 3 

x min {U w , 14.5} + 0.027 (H s - 0.4) 



1.25-0.25 



0.5 

Tc 



1.3 



(27) 



where f c is in Hertz, [7 10 is in meters per second, and H s 
is in meters. 

This relationship is well reproduced in hindcasts using 
C cu = —0.4 and s u = 1, while the BAJ source terms give a 
nearly constant value of 777,3 when H s varies and U\q is fixed 
(Ardhuin et al. 2009b). Here we also consider the fourth 
moment 7714 which, for linear waves, is proportional to a 
surface mean square slope filtered at the frequency f c . Fig- 
ure 8 shows that for any given wind speed mssc increases 
with the wave height (Gourrion et al. 2002), whereas this 
is not the case of 7774 in the BAJ parameterization, or, 
for very high winds, when C cu is too strong. In the case 
of BAJ, this is due to the (k/k r ) 2 part in the dissipation 
term (eq. 2), which plays a role similar to the cumula- 
tive term in our formulation. For C cu = — 1 and s u = 0, 
the cumulative effect gets too strong for wind speeds over 
10 m s _1 , in which case 7774 starts to decrease with in- 
creasing wave height, whereas for high winds and low (i.e. 
young) waves, the high frequency tail is too high and the 
mean square slope gets as large as 6%, which is unrealistic. 
It thus appears, that the high frequency tail, for s u = 0, 
responds too much to the wind, hence our use of s u = 1 in 
the TEST441 combination. The presence of a cumulative 
dissipation term allows a different balance in the spectral 
regions above the peak, where an equilibrium range with 
a spectrum proportional to / _4 develops (Long and Resio 
2007) , and in the high frequency tail were the spectrum de- 
cays like f~ 5 or possibly a little faster. The spectral level 
in the range 0.2 to 0.4 Hz was carefully compared against 
buoy data, and was found to be realistic. 

These interpretations of the model result assumes that 
the high frequency part of the spectrum can be simply con- 
verted to a wavenumber spectrum, using linear wave the- 
ory. This is not exactly the case as demonstrated by Ban- 
ner et al. (1989). Also, there is no consensus on the nature 
of the spectrum modelled with the energy balance equation 
but, since non-resonant nonlinearities are not represented, 
the modelled spectra are expected to be more related to 
Lagrangian buoy measurements, rather than Eulerian mea- 
surements. This matter is left for further studies, together 
with a detailed interpretation of altimeter radar cross sec- 



tions. Although it covers much less data, the analysis of 7714 
obtained from buoy heave spectra produces results similar 
to figure 8. 
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Fig. 8. Variation of the surface mean square slope es- 
timated as either 0.64/ero using the C-band altimeter on 
board JASON-1, after the correction of a 1.2 dB bias in 
the JASON data, or by integration of modelled spectra 
from to 0.72 Hz, with either the C cu = —0.4 and s u = 1 
parameterization (TEST441) or the parameterization by 
BAJ. For modelled values a constant 0.011 is added to ac- 
count for the short waves that contribute to the satellite 
signal and that are not resolved in the model. This satu- 
rated high frequency tail is consistent with the observations 
of Vandemark et al. (2004). The original 1 Hz data from 
JASON is subsampled at 0.5 Hz and averaged over 10 s, 
namely 58 km along the satellite track. The same averag- 
ing is applied to the wave model result, giving the 393382 
observations reported here, for the first half year of 2007. 



4. Verification 

In order to provide simplified measures of the difference 
between model time series A mo d and observations A" Q b s we 
use the following definitions for the normalized root mean 
square error (NRMSE), 



NRMSE(A) 



x„ 



J2 X obs 



(28) 
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the normalized bias, 



NB(X) = 



obs 



x 



mod 



obs 



(29) 



and Pearson's linear correlation coefficient, 



, y\ _ S (^obs — ^obs) (^mod — ^mod) , , 

V S (^obs ~ X b s ) (^mod — -^mod) 

where the overbar denotes the arithmetic average. 

The normalisation of the errors allows a quantitative 
comparison between widely different sea state regimes. Be- 
cause previous studies have often used (non-normalized) 
RMSE we also provide RMSE values. In addition to the 
coastal fetch-limited case of SHOWEX, presented above, 
the paramctcrizations are calibrated on at the global scale 
and validated in two other cases. 

a. Global scale results 

We present here results for the entire year 2007, using 
a stand-alone 0.5° resolution grid, covering the globe from 
80° south to 80° north. The model has actually been ad- 
justed to perform well over this data set, but the very large 
number of observations (over 2 million altimeter colloca- 
tion points) makes the model robust, and an independent 
validation on 2008 gives identical results. The interested 
reader may also look at the monthly reports for the SHOM 
model (e.g. Bidlot 2008), generated as part of the model 
verification project of the IOC-WMO Joint Commission 
on Oceanography and marine Meteorology (JCOMM), in 
which the TEST441 parameterization (C cu = —0.4 and 
s u = 1) is used, except for the Mediterranean where, the 
TEST405 has been preferred for its superior performance 
for younger seas. These SHOM models are ran in a combi- 
nation of two-way nested grids (Tolman 2007) . The monthly 
JCOMM reports include both analysis and forecasts, but, 
since they are produced in a routine setting, many SHOM 
calculations from December 2008 to June 2009 have been 
affected by wind file transfer problems. 

Comparing model results for H s to well-calibrated (Qucf- 
fculou and Croize-Fillon 2008) altimeter-derived measure- 
ments provides a good verification of the model perfor- 
mance in a number of different wave climates. Figure (9) 
shows that, as expected, the important positive bias in the 
swell-dominated regions when using the BAJ parameter- 
ization, has been largely removed. This is essentially the 
signature of the specific swell dissipation that is parameter- 
ized in S OU f The largest bias pattern now appears in the 
southern ocean, reaching 30 cm in the Southern Atlantic. 
Although this bias is small compared to the local averaged 
wave height, it is rather strange when the model errors are 
plotted as a function of wave height in figure (10). Why 



would the model overestimate the Southern ocean waves 
but understimate the very large waves? 

The structure of the large bias, also seen in model re- 
sults with BAJ, is reminiscent of the observed pattern in 
iceberg distribution observed by Tournadre et al. (2008). 
These observed iceberg distributions are enough to give a 
cross-section for incoming waves of the order of 1 to 10% for 
a 250 km propagation length. Taking icebergs into account 
could actually reverse the sign of the bias. This matter will 
be investigated elsewhere. 

Also noticeable is a significant negative bias in the equa- 
torial south Pacific, amplified from the same bias obtained 
with the BAJ parameterization. It is possible that the 
masking of subgrid islands (Tolman 2003) introduces a 
bias by neglecting shoreline reflections. This model de- 
fect could be exacerbated in this region by the very large 
ratio of shoreline length to sea area. This will also require 
further investigation. Finally, the negative biases for H s 
on mid-latitude east coasts are reduced but still persist. 
It is well known that these areas are also characterized 
by strong boundary currents (Gulf Stream, Kuroshio, Ag- 
ulhas ...) with warm waters that is generally conducive 
to wave amplification and faster wind-wave growth (e.g. 
Vandemark et al. 2001). Neither effect is included in the 
present calculation because the accuracy of both modelled 
surface currents and air-sea stability paramctcrizations are 
likely to be insufficient (Collard et al. 2008; Ardhuin et al. 
2007). 

The reduction of systematic biases clearly contributes 
to the reduction of r.m.s. errors, as evident in the equa- 
torial east Pacific (Fig. 11). However, the new parame- 
terization also brings a considerable reduction of scatter, 
with reduced errors even where biases are minimal, such as 
the trade winds area south of Hawaii, where the NRMSE 
for H s can be as low as 5%. When areas within 400 km 
from continents are excluded, because the global model 
resolution may be inadequate, significant errors (> 12.5%, 
in yellow to red) remain in the northern Indian ocean, on 
the North American and Asian east coasts, the Southern 
Atlantic. The paramctcrizations TEST405, TEST437 and 
TEST441 produce smaller errors on average than BAJ. It 
is likely that the model benefits from the absence of swell 
influences on wind seas: swell in BAJ typically leads to a 
reduced dissipation and stronger wind wave growth. As 
models are adjusted to average sea state conditions, this 
adjustement leads to a reduced wind sea growth on east 
coasts where there is generally less swell. 

Although much more sparse than the altimeter data, 
the in situ measurements collected and exchanged as part 
as the JCOMM wave model verification is very useful for 
constraining other aspects of the sea state. This is illus- 
trated here with mean periods T m o2 for data provided by 
the U.K. and French meteorological services, and peak pe- 
riods T p for all other sources. It is worth noting that the 
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Fig. 9. Bias for the year 2007 in centimeters. The global 0.5 WWATCH model is compared to altimeters JASON, EN- 
VIS AT and GFO following the method of Rascle et al. (2008). The top panel is the result with the BAJ parameterization, 
and the bottom panel is the result with the C cu = —0.4 and s u — 1 (TEST441) parameterization. 
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Normalized RMSE (%) 

Fig. 11. Normalized RMSE for the significant wave height over the year 2007, in percents. The global 0.5° resolution 
WWATCH model is compared to altimeters JASON, ENVISAT and GFO following the method of Rascle et al. (2008). 
The top panel is the result with the BAJ parameterization, and the bottom panel is the result with the C cu = —0.4 and 
s u = 1 (TEST441) parameterization. 
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Fig. 10. Wave model errors as a function of H s . All model 
parameterizations are used in a global WWATCH model 
settings using a 0.5° resolution. The model output at 3h 
intervals is compared to JASON, ENVISAT and GFO fol- 
lowing the method of Rascle et al. (2008). Namely, the 
altimeter 1 Hz Ku band estimates of H s are averaged over 
1°. After this averaging, the total number of observations 
is 2044545. The altimeter estimates are not expected to be 
valid for H s larger than about 12 m, due to the low signal 
level and the fact that the waveform used to estimate H s 
is not long enough in these cases. 



dissipation was tested but it tended to deteriorate the re- 
sults on other parameters. On European coasts, despite a 
stronger bias, the errors on T m o2 are particularly reduced. 
Again, this reduction of the model scatter can be largely 
attributed to the decoupling of swell from windsca growth. 

Table 1. Model accuracy for measured wave parameters 
over the oceans in 2007. mss data from JASON 1 corre- 
sponds to January to July 2007 (393382 co- located points). 
Unless otherwise specified by the number in parenthesis, 
the cut-off frequency is take to be 0.4 Hz, C stands for 
C-band. The normalized bias (NB) is defined as the bias 
divided by the r.m.s. observed value, while the scatter in- 
dex (SI) is defined as the r.m.s. difference between modeled 
and observed values, after correction for the bias, normal- 
ized by the r.m.s. observed value, and r is Pearson's corre- 
lation coefficient. These global averages are area- weighted, 
and the SI and NRMSE are the area-weighted averages of 
the local SI and NRMSE. 
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errors on H s for in situ platforms arc comparable to the 
errors against altimeter data. 

With the BAJ parameterizations, the largest errors in 
the model results are the large biases on peak periods on 
the U.S. West coast (Fig. 12), by 1.2 to 1.8 s for most loca- 
tions, and the understimation of peak periods on the U.S. 
East Coast. However, peak and mean periods off the Euro- 
pean coasts were generally very well predicted. When using 
the TEST441 parameterization, the explicit swell dissipa- 
tion reduces the bias on periods on the U.S. West coast, 
but the problem is not completely solved, with residual bi- 
ases of 0.2 to 0.4 s. This is consistent with the validation 
using satellite SAR data (Fig. 1), that showed a tendency 
to underpredict steep swells near the storms and overpre- 
dict them in the far field. A simple increase of the swell 



The general performance of the parameterizations is 
synthetized in table 1. It is interesting to note that the 
parameterization TEST405 that uses a diagnostic tail for 
2.5 times the mean frequency gives good results in terms 
of scatter and bias even for parameters related to short 
waves (TO3, mi). This use of diagnostic tail is thus a good 
pragmatic alternative to the more costly explicit resolu- 
tion of shorter waves, which requires a smaller adapata- 
tive timestep, and more complex parameterizations. The 
diagnostic tail generally mimics the effect of both the cu- 
mulative and sheltering effects. Yet, the parameterization 
TEST441 demonstrates that it is possible to obtain slightly 
better results with a free tail. The normalized biases indi- 
cated for the mean square slopes are only relative because 
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Fig. 12. Statistics for the year 2007 based on the JCOMM verification data base (Bidlot et al. 2007b). Bais (a,d) and 
NRMSE (b,c) for T p or T m o2 at in situ locations using the BAJ or the proposed TEST441 (T p is shown at all buoys except 
U.K. and French buoys for which T m o2 is shown). The different symbols are only used to help distinguish the various 
colors and do not carry extra information. 
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of the approximate calibration of the radar cross section. 
They show that the BAJ parameterization (Bidlot et al. 
2005), and to a lesser extent the use of a f~ 5 tail, produce 
energy levels that are relatively lower at high frequency. 

b. Lake Michigan 

At the global scale, the sea state is never very young, 
and it is desirable to also verify the robustness of the pa- 
rameterization in conditions that are more representative 
of the coastal ocean. We thus follow the analysis of wave 
model performances by Rogers and Wang (2007), here- 
inafter RW2007, and give results for the Lake Michigan, 
representative of relatively young waves. The model was 
applied with parameterizations BAJ, TC, TEST437 and 
TEST441 over the same time frame as investigated by 
RW2007, September 1 to November 14, 2002. The model 
setting and forcing fields are identical to the one defined 
by (Rogers et al. 2003), with a 2 km resolution grid, a 10° 
directional resolution, and a wind field defined from in situ 
observations. The results at the position of National Data 
Buoys Center's buoy 45007 are compared to corresponding 
measurements. 

Using the directional validation method proposed by 
these authors, the TC parameterization undcrstimates the 
directional spread ag by 1.2 to 1.6° in the range 0.8 to 2.0 
f p , and more at higher frequencies. The understimation 
with BAJ is about half, and the TEST441 and TEST437 
overpredict the directional spread by about 2.3 to 5.9° in 
the range 0.8 to 2.0/ p , and less so for higher frequencies. 
It thus appears that the broadening introduced to fit the 
SHOWEX 1999 observations is not optimal for other situ- 
ations. A similar positive bias on directional spread is also 
found in global hindcasts. 

Further results are presented in Table 2 and Figure 13. 
The top panel of Figure 13 compare the summed values 
of co-located model and observed spectral densities for the 
duration of the simulation. This presentation provides fre- 
quency distribution of bias of the various models, while also 
indicating the relative contribution of each frequency to the 
wave climate for this region and time period. The lower 
panel shows the correlation coefficient r for the equivalent 
significant waveheights computed for multiple frequency 
bands. This is presented in terms of f j fp (bin width=0.1), 
with f p being calculated as the stabilized "synthetic peak 
frequency" of the corresponding buoy spectrum, as defined 
in RW2007. 

The most noticeable outcome of these comparisons is 
the relatively poor performance of the TC parameteriza- 
tions. Taken in context with other TC results presented 
herein and prior undocumented application of the model in 
the Great Lakes with model wind fields, this suggests that 
these parameterizations have some undesirable dependence 
on scale, with the parameters adjusted by Tolman (2002b) 
being most optimal for ocean-scale applications. 



The KHH parameterization perform well in this simple 
wave climate, consistent with prior published applications 
with the SWAN model, (Rogers et al. 2003) and RW2007, 
without the difficulties discussed above in mixed sea-swell 
conditions. 

The BAJ, TEST437 and TEST441 models also perform 
well here. Taken together with the global comparisons 
above, we observe no apparent dependency of model skill 
on the scale of the application with these three physics. In 
the bias comparison, Fig. 13 top panel, the BAJ model 
is nearly identical to the KHH model. Similarly, the two 
new models are also very close. Although TEST441 and 
TEST437 produce a minor underestimation of H s (Table 
2), they give slightly better correlations with observed wave 
heights and mean periods. 

Table 2. Model-data comparison at NDBC buoy 45007 
(Lake Michigan) for four frequency bands. Statistics are 
given for equivalent significant wave heights in the bands 
0.5/p < f < 0.8/p (band 1), 0.8/ p < f < 1.2 f p (band 2), 
1.2/ p < / < 2f p (band 3) 2f p < f < 3/ p (band 4). The 
KHH n = 2 run corresponds to the dissipation parameter- 
ization defined by (Rogers et al. 2003) based on (Komen 
et al. 1984) and applied in the WW3 code. For the peak 
frequency band, statistics are also given for the directional 
spread ag. 
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It thus appears that for such young seas, the directional 
spreading of the parameterization could be improved, but 
the energy content of various frequency bands, and as a 
result the mean period, are reproduced with less scatter 
than with previous parameterizations. 
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Fig. 13. Model-data comparison at NDBC buoy 45007 
(Lake Michigan). Top panel: Comparison of summed spec- 
tral density versus frequency for the duration of the sim- 
ulation. Lower panel: Correlation coefficients versus nor- 
malized frequency (see text for explanation). 



c. Hurricane Ivan 

Although the global hindcast does contain quite a few 
extreme events, with significant wave heights up to 17 m, 
these were obtained with a relatively coarse wave model 
grid and wind forcing (0.5° resolution and 6-h timestep) 
that is insufficient to resolve small storms such as tropical 
cyclones (Tolman and Alves 2005). Hurricane waves do 
share many similarities with more usual sea states (Young 
2005), but the high winds and their rapid rotation are par- 
ticularly challenging for numerical wave models. It is thus 
necessary to verify that the new source functions perform 
adequately under extreme wind conditions. A simulation 
of Hurricane Ivan (Gulf of Mexico, September 2004) is cho- 
sen for this purpose because it was extensively measured 
(e.g. Wang et al. 2005) and hindcastcd. 

Winds for this simulation are based on gridded surface 
wind analyses created by NOAA's Hurricane Research Di- 
vision (HRD). These analyses are at three hour intervals, 
which for a small, fast-moving weather system is tempo- 
rally too coarse to provide directly to the wave model. 
Therefore, as an intermediate step, fields are reprocessed 
to 30 minute intervals, with the storm position updated at 
each interval (thus, semi-Lagrangian interpolation). The 
wind speeds are reduced by factor 1/1.11 to convert from 
maximum sustained gust to hourly mean. The HRD winds 
do not cover the entire computational domain. For areas 
falling outside the domain, the nearest NDBC buoy wind 
observation is used. This produces some non-physical spa- 
tial discontinuities in the wind field, but these are smoothed 



in the wave model integration, and in any event, only affect 
weaker wind seas far from the storm center. 

Bathymetry is taken from the Naval Research Labora- 
tory's 2 minute resolution database, DBDB2, coarsened to 
the model grid resolution (0.1 deg). The directional reso- 
lution is 10°, and the frequency range is limited to 0.0418- 
0.4117 Hz. The model was applied from September 13 to 
September 16 2004. Model results are illustrated by figure 
14. The models were validated at all the buoys in the gulf 
of Mexico. Results from buoy 42040, where waves where 
largest, are shown here. 
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Fig. 14. Time series of model and buoy significant wave 
height (top) and partial wave height (bottom, for / < 0.06 
Hz), during the passage of Hurrican Ivan, at buoy 42040. 
The 10-m discus buoy was capized by waves and could not 
record waves after 16 September. 

Model runs with paramctcrizations BAJ, TEST437 and 
TEST441 give very similar results: close to the observa- 
tions, except for the highest waves (H s > 13 m at buoy 
42040) where TEST437 and TEST441 give slightly smaller 
values. Results with the TC parameterization are generally 
lower in terms of H s than all these parameterizations that 
share a Janssen-type input. It appears that the new source 
term perform similarly to BAJ and are able to reproduce 
such young waves and severe sea states. 

Because the wind forcing enters the wave model through 
a wind stress parameterized in a way that may not apply 
to such conditions, it is worthwhile to re-examine some 
choices made above. In particular the surface roughness 
was allowed to exceed 0.002 in TEST441b, which resulted 
in better estimates of H s . Lifting this constraint shows 
that, for these very high winds, the wind sheltering ef- 
fect plays a similar role to the limitation of the roughness 
to zo.max, with the difference that it tends to narrow the 
wind input spectrum (Fig. 7). This narrower wind input 
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has a limited effect of the wind stress and H s , but is has a 
noticeable effect on the spectral shape. This is illustrated 
by the low frequency energy that appears to be strongly 
overestimated before the peak of the storm for TEST441b. 
In general the new parameterizations provide results that 
are as reasonable as those of previous parameterizations, 
given the uncertainty of the wind forcing. 

5. Conclusions 

A set of parameterization for the dissipation source 
terms of the wave energy balance equation have been pro- 
posed, based on known properties of swell dissipation and 
wave breaking statistics. This dissipation includes an ex- 
plicit nonlinear swell dissipation and a wave breaking pa- 
rameterization that contains a cumulative term, represent- 
ing the dissipation of short waves by longer breakers, and 
different dissipation rates for different directions. These 
dissipation parameterizations have been combined with a 
modified form of the wind input proposed by Janssen (1991), 
in which the questionable gustincss parameter z a has been 
reduced, and the general shape of the wind input has been 
significantly modified. The resulting source term balance is 
thus markedly different from the previous proposed forms, 
with a near-balance for very old seas between the air-sea 
friction term, that dissipates swell, and the nonlinear en- 
ergy flux to low frequencies. Also, the wind input is con- 
centrated in a narrower range of frequencies. 

For younger seas the wind input is relatively weaker 
than given by Janssen (1991) but stronger than given by 
Tolman and Chalikov (1996) (TC). However, the dissipa- 
tion at the peak is generally stronger because it is essen- 
tially based on a local steepness and these dominant waves 
are the steepest in the sea state. As a result the short 
fetch growth is relatively weaker than with the source term 
combination proposed by Bidlot ct al. (2007a) (BAJ). The 
choice of parameters tested here tend to produce broader 
directional spectra than observed in the Lake Michigan and 
global hindcasts, and slanting fetch directions that are too 
oblique relative to the wind (Fig. 2). In this respect the 
new source terms are intermediate between BAJ and TC. 

Another likely defect comes from the definition of the 
saturation level used to define the breaking-induced dissi- 
pation. Here, as in the work by van der Westhuysen et al. 
(2007), the saturation is local in frequency space, whereas 
wave breaking is naturally expected to have a relatively 
broad impact due to its localization in space and time 
(Hasselmann 1974). This is expected to produce an overes- 
timation of the energy just below the peak, and an under- 
stimation at the peak of the saturation spectrum. These 
effects likely contribute to the persistent overestimation of 
low frequency energy in the model. 

In spite of these defects, the new parameterization pro- 
duces robust results and clearly outperform the Bidlot et al. 



(2007a) parameterization in global hindcasts, whether one 
considers dominant wave parameters, H Sl T m o2 and T p or 
parameters sensitive to the high frequency content, such as 
the surface Stokes U ss drift or the mean square slope. At 
global scales, errors on iJ s , T p and U ss are - on average - re- 
duced by 15, 25 and 50% relative to those obtained with the 
parameterization by Bidlot ct al. (2007a). All global hind- 
casts results, for the year 2007 at least, are available for fur- 
ther analysis at ftp://ftp.ifremer.fr/ifrcmer/cersat/products/griddcd, 

Another important aspect was the validation at regional 
to global scales. We note that the TC paper does include 
verification with steady-state, fetch-limited growth curves. 
Though such verification is a useful step, the outcome of 
the Lake Michigan hindcast suggests that such verification 
does in no way anticipate skill in real sub-rcgional-scale 
applications. One of the parameterizations proposed here 
(TEST441) also gives slightly lower performance for young 
seas, which is not obvious in the case of lake Michigan, 
but was revealed by hindcasts of Mediterranean waves (not 
shown) . 

Because our intention was only to demonstrate the ca- 
pability of new dissipation parameterizations and the re- 
sulting source term balances, we have not fully adjusted 
the 18 parameters that define the deep water parameter- 
izations, compared to about 9 with Bidlot ct al. (2007a). 
The results presented here are thus preliminary in terms of 
model performance, which is why the parameterizations are 
still given temporary names like TEST441. As illustrated 
by the hurrican Ivan hindcast, some parameters, such as 
zo,max, is probably unnecessary: in that particular case the 
removal of Zo.max improved the results, but for global scale 
results it had no impact at all (not shown). 

Because 5 of the extra parameters define the air-sea 
friction term that produces swell dissipation, and 2 de- 
fine the cumulative breaking term, it is feasible to define 
a systematic adjustement procedure that should produce 
further improvements by separately adjusting swell, wind 
sea peak, and high frequency properties. In particular the 
directional distribution may be improved by making the 
dissipation term more isotropic (i.e. taking 84 > 0.3) or 
modifying the definition of the saturation parameter B' in 
equation (12). In part II we shall further investigate the 
response of the wave field to varying currents, from global 
scales to regional tidal currents. It is particularly expected 
that wave steepening will produce much more dissipation 
due to breaking, as envisaged by Phillips (1984). 

Obviously, it is well known that the Discrete Interac- 
tion Approximation used here to compute the non-linear 
interactions is the source of large errors, and further calcu- 
lations, will be performed using a more accurate estimation 
of these interactions in part III. 
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APPENDIX A 

Parameter settings 

All parameters defining the dissipation source function 
and their numerical values are listed in table 3 for the wind- 
wave interaction term S a tm and table 4 for the wave-ocean 
interaction term S oc . We also recall that the nonlinear cou- 
pling coefficient (variable NLPROP in WWATCH) is set to 
2.78 x 10 7 in all cases, except for the two parameterizations 
mostly used here, with C„i = 2.5 x 10 7 in TEST437 and 
TEST441. Although the best performance for most param- 
eters is obtained with the TEST441 settings, its underes- 
timation of extreme sea states may be a problem in some 
applications for which the TEST437 may be preferred. A 
full tuning of the model has not been tried yet and it is 
possible that a simple adjustment of /3 ma x C cu , r cu and s u 
may produce even better results. Finally, these parame- 
ters have been mostly adjusted for deep water conditions 
using ECMWF winds. Using other sources of winds for 
large scale applications may require a retuning of the wind 
source function, which can be best performed by a read- 
justment of /3 max . 



REFERENCES 

Alves, J. H. G. M. and M. L. Banner, 2003: Performance of 
a saturation-based dissipation-rate source term in mod- 
eling the fetch-limited evolution of wind waves. J. Phys. 
Oceanogr., 33, 1274-1298. 

Alves, J. H. G. M., M. L. Banner, and I. R. Young, 2003: 
Revisiting the Pierson-Moskowitz asymptotic limits for 
fully developed wind waves. J. Phys. Oceanogr., 33, 
1301-1323. 

Ardhuin, F., B. Chapron, and F. Collard, 2009a: Obser- 
vation of swell dissipation across oceans. Geophys. Res. 
Lett, 36, doi:10.1029/2008GL037030, L06607. 



Ardhuin, F., F. Collard, B. Chapron, P. Queffeulou, J.-F. 
Filipot, and M. Hamon, 2008a: Spectral wave dissipation 
based on observations: a global validation. Proceedings 
of Chinese- German Joint Symposium on Hydraulics and 
Ocean Engineering, Darmstadt, Germany, 391-400. 

Ardhuin, F., T. H. C. Herbers, K. P. Watts, G. P. van Vled- 
der, R. Jensen, and H. Graber, 2007: Swell and slanting 
fetch effects on wind wave growth. J. Phys. Oceanogr., 
37, doi:10.1175/JPO3039.1, 908-931. 

Ardhuin, F. and A. D. Jenkins, 2005: On the effect of 
wind and turbulence on ocean swell. Proceedings of the 
15th International Polar and Offshore Engineering Con- 
ference, June 19-24, Seoul, South Korea, ISOPE, vol- 
ume III, 429-434. 

- 2006: On the interaction of surface waves and upper 
ocean turbulence. J. Phys. Oceanogr., 36, 551-557. 

Ardhuin, F. and A. Le Boyer, 2006: Numerical modelling 
of sea states: validation of spectral shapes (in French). 
Navigation, 54, 55-71. 

Ardhuin, F., L. Marie, N. Rasclc, P. Forget, and A. Roland, 
2009b: Observation and estimation of Lagrangian, 
Stokes and Eulerian currents induced by wind and waves 
at the sea surface. J. Phys. Oceanogr., 39, 2820-2838. 

Ardhuin, F., N. Rascle, and K. A. Bclibassakis, 2008b: 
Explicit wave-averaged primitive equations using a 
generalized Lagrangian mean. Ocean Modelling, 20, 
doi:10.1016/j.ocemod.2007.07.001, 35-60. 

Babanin, A., K. Tsagareli, I. Young, and D. Walker, 2007: 
Implementation of new experimental input/dissipation 
terms for modelling spectral evolution of wind waves. 
Proceedings of the 1 0th International workshop on wave 
hindcasting and forecasting, Oahu, Hawaii. 

Babanin, A., I. Young, and M. Banner, 2001: Breaking 
probabilities for dominant surface waves on water of fi- 
nite depth. J. Geophys. Res., 106, 11659-11676. 

Babanin, A. V. and A. J. van der Westhuysen, 2008: 
Physics of saturation-based dissipation functions pro- 
posed for wave forecast models. J. Phys. Oceanogr., 38, 
1831-1841. 

Babanin, A. V. and I. R. Young, 2005: Two-phase 
behaviour of the spectral dissipation of wind waves. 
Proceedings of the 5th International Symposium Ocean 
Wave Measurement and Analysis, Madrid, june 2005, 
ASCE, paper number 51. 

Banner, M. L., A. V. Babanin, and I. R. Young, 2000: 
Breaking probability for dominant waves on the sea sur- 
face. J. Phys. Oceanogr., 30, 3145-3160. 



21 



Banner, M. L., J. R. Gemmrich, and D. M. Farmer, 2002: 
Multiscale measurement of ocean wave breaking proba- 
bility. J. Phys. Oceanogr., 32, 3364-3374. 

Banner, M. L., I. S. F. Jones, and J. C. Trinder, 1989: 
Wavenumber spectra of short gravity waves. J. Fluid 
Mech., 198, 321-344. 

Banner, M. L. and R. P. Morison, 2006: On modeling 
spectral dissipation due to wave breaking for ocean wind 
waves. Proceedings of the 9th International workshop on 
wave hindcasting and forecasting, Victoria, Canada. 

Banner, M. L. and W. L. Peirson, 2007: Wave breaking 
onset and strength for two-dimensional deep-water wave 
groups. J. Fluid Mech., 585, 93-115. 

Banner, M. L. and I. R. Young, 1994: Modeling spectral 
dissipation in the evolution of wind waves, part I: assess- 
ment of existing model performance. J. Phys. Oceanogr., 
24, 1550-1570. 

Barber, N. F., 1949: Behaviour of waves on tidal streams. 
Proc. Roy. Soc. Lond. A, 198, 81-93. 

Bidlot, J., S. Abdalla, and P. Janssen, 2005: A revised for- 
mulation for ocean wave dissipation in CY25R1. Techni- 
cal Report Memorandum R60.9/JB/0516, Research De- 
partment, ECMWF, Reading, U. K. 

Bidlot, J., P. Janssen, and S. Abdalla, 2007a: A revised for- 
mulation of ocean wave dissipation and its model impact. 
Technical Report Memorandum 509, ECMWF, Reading, 
U. K. 

Bidlot, J.-R., 2008: Intercomparison of operational wave 
forecasting systems against buoys: data from ECMWF, 
MetOffice, FNMOC, NCEP, DWD, BoM, SHOM and 
JMA, September 2008 to November 2008. Technical 
report, Joint WMO-IOC Technical Commission for 
Oceanography and Marine Meteorology, available from 
http: / /preview. tinyurl.com/7bz6jj. 

Bidlot, J.-R., J.-G. Li, P. Wittmann, M. Fauchon, H. Chen, 
J.-M. Lefevre, T. Bruns, D. Greenslade, F. Ardhuin, 
N. Kohno, S. Park, and M. Gomez, 2007b: Inter- 
comparison of operational wave forecasting systems. Pro- 
ceedings, 10th Int. WOrkshop of Wave Hindcasting and 
Forecasting, Hawaii. 

Cavaleri, L., 2006: Wave modeling where to go in the fu- 
ture. Bull. Amer. Meterol. Soc, 87, 207-214. 

Chalikov, D. V. and M. Y. Belevich, 1993: Onc- 
dimcnsional theory of the wave boundary layer. 
Boundary-Layer Meteorol, 63, 65-96. 



Chen, G. and S. E. Belcher, 2000: Effects of long waves 
on wind-generated waves. J. Phys. Oceanogr., 30, 2246- 
2256. 

Collard, F., F. Ardhuin, and B. Chapron, 2009: Monitor- 
ing and analysis of ocean swell fields using a spacebornc 
SAR: a new method for routine observations. J. Geophys. 
Res., 114, C07023. 

Collard, F., A. Mouche, B. Chapron, C. Danilo, 
and J. Johannessen, 2008: Routine high resolu- 
tion observation of selected major surface cur- 
rents from space. Proceedings of SEASAR 2008, 
SP-656, ESA, ESA - ESRIN, Frascati, Italy, 
http: / / earth. esa.int /workshops / seasar2008 /participants /287/ pres. 

Dalrymplc, R. A., 1974: A finite amplitude wave on a linear 
shear current. J. Geophys. Res., 79, 4498-4504. 

Darbyshire, J., 1958: The generation of waves by wind. 
Phil. Trans. Roy. Soc. London A, 215, 299-428. 

Dobson, F., W. Perric, and B. Toulany, 1989: On the 
deep water fetch laws for wind-generated surface gravity 
waves. Atmosphere Ocean, 27, 210-236. 

Donelan, M. A., 1998: Air-water exchange processes. Phys- 
ical Processes in Lakes and Oceans, J. Imbcrgcr, ed., 
American Geophysical Union, Washington, D.C., pages 
18-36, iSBN 0-87590-268-5. 

Dore, B. D., 1978: Some effects of the air- water interface 
on gravity waves. Geophys. Astrophys. Fluid. Dyn., 10, 
215-230. 

Filipot, J.-F., F. Ardhuin, and A. Babanin, 2010: A unified 
deep-to-shallow-water spectral wave breaking dissipation 
formulation. Part 1. Breaking probability. J. Geophys. 
Res., 115, in press. 

Gclci, R., H. Cazale, and J. Vassal, 1957: Prevision de 
la houlc. La methode des densites spectroangulaircs. 
Bulletin d 'information du Comite d'Oceanographie et 
d'Etude des Cotes, 9, 416-435. 

Gemmrich, J. R., M. L. Banner, and C. Garrett, 2008: 
Spectrally resolved energy dissipation rate and momen- 
tum flux of breaking waves. J. Phys. Oceanogr., 38, 
1296-1312. 

Gourrion, J., D. Vandemark, S. Bailey, and B. Chapron, 
2002: Investigation of C-band altimeter cross section de- 
pendence on wind speed and sea state. Can. J. Remote 
Sensing, 28, 484-489. 

Grant, W. D. and O. S. Madsen, 1979: Combined wave and 
current interaction with a rough bottom. J. Geophys. 
Res., 84, 1797-1808. 



22 



Hargreaves, J. C. and J. D. Annan, 2000: Comments on 
improvement of the short-fetch behavior in the wave 
ocean model (WAM)-. J. Atmos. Ocean Technol., 18, 
711-715. 

Harris, D. L., 1966: The wave-driven wind. J. Atmos. Sci., 
23, 688-693. 

Hassclmann, K., 1971: On the mass and momentum trans- 
fer between short gravity waves and larger-scale motions. 
J. Fluid Mech., 4, 189-205. 

— 1974: On the spectral dissipation of ocean waves due to 
white capping. Boundary-Layer Meteoroi, 6, 107-127. 

Hassclmann, S., K. Hassclmann, J. Allender, and T. Bar- 
nett, 1985: Computation and paramctcrizations of the 
nonlinear energy transfer in a gravity-wave spectrum. 
Part II: Paramctcrizations of the nonlinear energy trans- 
fer for application in wave models. J. Phys. Oceanogr., 
15, 1378-1391. 

Hogstrom, U., A. Smcdman, E. Sahlee, H. Pettersson and 

F. Zhang, 2009: H. Pettersson, and F. Zhang, 2009: The 
atmospheric boundary layer during swell: A held study 
and interpretation of the turbulent kinetic energy budget 
for high wave ages. J. Atmos. Sci., 66, 2764-2779. 

Jansscn, P. A. E. M., 1991: Quasi-linear theory of wind 
wave generation applied to wave forecasting. J. Phys. 
Oceanogr., 21, 1631-1642, see comments by D. Chalikov, 
J. Phys. Oceanogr. 1993, vol. 23 pp. 1597-1600. 

Janssen, P. A. E. M., K. Hassclmann, S. Hasselmann, and 

G. J. Komen, 1994: Parameterization of source terms 
and the energy balance in a growing wind sea. Dynamics 
and modelling of ocean waves, G. J. Komen et al., ed., 
Cambridge University Press, pages 215-238. 

Jansscn, P. A. E. M., O. Saetra, C. Wettre, and H. Hers- 
bach, 2004: Impact of the sea state on the atmosphere 
and ocean. Annates Hydrographiques , 6e serie, vol. 3, 
3-1-3-23. 

Jensen, B. L., B. M. Sumer, and J. Freds0e, 1989: Turbu- 
lent oscillatory boundary layers at high Reynolds num- 
bers. J. Fluid Mech., 206, 265-297. 

Komen, G. J., L. Cavalcri, M. Donelan, K. Hasselmann, 
S. Hassclmann, and P. A. E. M. Jansscn, 1994: Dynam- 
ics and modelling of ocean waves. Cambridge University 
Press, Cambridge, 554 pp. 

Komen, G. J., K. Hasselmann, and S. Hasselmann, 1984: 
On the existence of a fully developed windsea spectrum. 
J. Phys. Oceanogr., 14, 1271-1285. 



Lefevre, J.-M., S. E. §tefanescu, and V. Makin, 2004: Im- 
plementation of new source terms in a third generation 
wave model. Preprints of the 3th International workshop 
on wave hindcasting and forecasting, Montreal, Quebec, 
19-22 May, Environment Canada. 

Long, C. E. and D. T. Resio, 2007: Wind wave spectral ob- 
servations in Currituck Sound, North Carolina. J. Geo- 
phys. Res., 112, doi:10.1029/2006JC003835, C05001. 

Longuet-Higgins, M. S. and J. S. Turner, 1974: An 'en- 
training plume' model of a spilling breaker. J. Fluid 
Mech., 63, 1-20. 

Magne, R., K. Belibassakis, T. H. C. Herbers, F. Ardhuin, 
W. C. O'Reilly, and V. Rey, 2007: Evolution of surface 
gravity waves over a submarine canyon. J. Geophys. Res., 
112, doi:10.1029/2005JC003035, C01002. 

Manasseh, R., A. V. Babanin, C. Forbes, K. Rickards, 
I. Bobevski, and A. Ooi, 2006: Passive acoustic determi- 
nation of wave-breaking events and their severity across 
the spectrum. J. Atmos. Ocean. Tech., 23(4), 599-618. 

Makin, V. K. and M. Stam, 2003: New drag formulation 
in NEDWAM. Technical Report 250, Koninklijk Ned- 
crlands Meteorologisch Instituut, Dc Bilt, The Nether- 
lands. 

Melville, W. K, F. Vcrron, and C. J. White, 2002: The 
velocity field under breaking waves: coherent structures 
and turbulence. J. Fluid Mech., 454, 203-233. 

Munk, W. H. and M. A. Traylor, 1947: Refraction of 
ocean waves: a process linking underwater topography 
to beach erosion. Journal of Geology, LV, 1-26. 

Peirson, W. L. and M. L. Banner, 2003: Aqueous surface 
layer flows induced by microscalc breaking wind waves. 
J. Fluid Mech., 479, 1-38. 

Phillips, O. M., 1958: The equilibrium range in the spec- 
trum of wind-gencrated waves. J. Fluid Mech., 4, 426- 
433. 

- 1963: On the attenuation of long gravity waves by short 
breaking waves. J. Fluid Mech., 16, 321-332. 

- 1984: On the response of short ocean wave components 
at a fixed wavenumber to ocean current variations. J. 
Phys. Oceanogr., 14, 1425-1433. 

- 1985: Spectral and statistical properties of the equi- 
librium range in wind-generated gravity waves. J. Fluid 
Mech., 156, 505-531. 

Polnikov, V. G. and V. Inoccntini, 2008: Comparative 
study of performance of wind wave model: Wavewatch 
modified by new source function. Engineering Applica- 
tions of Computational Fluid Mechanics , 2, 466-481. 



23 



Queffeulou, P. and D. Croize-Fillon, 2008: Tolman, H. L. and J.-H. G. M. Alves, 2005: Numerical 
Global altimeter SWH data set, version modeling of wind waves generated by tropical cyclones 
4, October 2008. Technical report, Ifremer. using moving grids. Ocean Modelling, 9, 305-323. 

ftp:/ /ftp. ifrcmcr.fr/ifremer/cersat/products/swath/altimeters/waves/ 
documentation/altimeter _wave_merge__4.0.pdf 



Rascle, N., F. Ardhuin, P. Queffeulou, and D. Croize- 
Fillon, 2008: A global wave parameter database for geo- 
physical applications, part 1: wave-current-turbulcncc 
interaction parameters for the open ocean based on tra- 
ditional paramctcrizations. Ocean Modelling, 25, 154- 
171, doi:10.1016/j.occmod.2008.07.006. 

Rogers, W. E., P. A. Hwang, and D. W. Wang, 2003: Inves- 
tigation of wave growth and decay in the SWAN model: 
Three regional-scale applications. J. Phys. Oceanogr., 
33, 366-389. 

Rogers, W. E. and D. W. C. Wang, 2007: Directional val- 
idation of wave predictions. J. Atmos. Ocean Technoi, 
24, 504-520. 

Ruessink, B. G., D. J. R. Walstra, and H. N. Southgatc, 
2003: Calibration and verification of a parametric wave 
model on barred beaches. Coastal Eng., 48, 139-149. 

Smcdman, A., U. Hogstrom, E. Sahlee, W. M. Drennan, 
K. K. Kahma, H. Pettersson and F. Zhang, 2009: Obser- 
vational Study of Marine Atmospheric Boundary Layer 
Characteristics during Swell. J. Atmos. Sci., 66, 2747- 
2763. 

Stansell, P. and C. MacFarlanc, 2002: Experimental inves- 
tigation of wave breaking criteria based on wave phase 
speeds. J. Phys. Oceanogr., 32, 1269-1283. 

Tolman, H. L., 1992: Effects of numerics on the physics in 
a third-generation wind- wave model. J. Phys. Oceanogr., 
22, 1095-1111. 

- 2002a: Limiters in third-generation wind wave models. 
Global Atmos. Ocean Syst., 8, 67-83. 

- 2002b: Validation of WAVEWATCH-IH version 1.15. 
Technical Report 213, NOAA/NWS/NCEP/MMAB. 

- 2003: Treatment of unresolved islands and ice in wind 
wave models. Ocean Modelling, 5, 219-231. 

- 2007: The 2007 release of WAVEWATCH EI. Proceed- 
ings, 10th Int. Workshop of Wave Hindcasting and Fore- 
casting, Hawaii. 

- 2008: A mosaic approach to wind wave modeling. 
Ocean Modelling, 25, doi:10.1016/j.ocemod.2008.06.005, 
35-47. 

- 2009: User manual and system documentation of 
WAVEWATCH-LU™ version 3.14. Technical Report 
276, NOAA/NWS/NCEP/MMAB. 



Tolman, H. L. and D. Chalikov, 1996: Source terms in a 
third-generation wind wave model. J. Phys. Oceanogr., 
26, 2497-2518. 



Tournadrc, J., K. Whitmcr, and F. Girard- Ardhuin, 
2008: Iceberg detection in open water by al- 
timeter waveform analysis. J. Geophys. Res., 113, 
doi: 10. 1029 /2007JC004587, C08040. 

Tsagareli, K., 2008: Numerical Investigation of Wind Input 
and Spectral Dissipation in Evolution of Wind Waves. 
Ph.D. thesis, University of Adelaide, Faculty of Engi- 
neering, Computer and Mathematical Sciences School of 
Civil and Environmental Engineering, Australia. 

Tulin, M. P. and M. Landrini, 2001: Breaking waves in 
the ocean and around ships. Proceedings of the 23rd 
ONR symposium on naval hydrodynamics, Vol de Reuil, 
France, Naval Studies Board, 713-745. 

van der Westhuyscn, A. J., M. Zijlema, and J. A. Bat- 
tjes, 2007: Saturation-based whitecapping dissipation in 
SWAN for deep and shallow water. Coastal Eng., 54, 
151-170. 

van der Westhuyscn, A. J., 2007: Advances in the spectral 
modelling of wind waves in the nearshore. Ph.D thesis, 
Delft University of Technology, 206 pp. 

van Vledder, G. P. and D. P. Hurdle, 2002: Performance of 
formulations for whitecapping in wave prediction models. 
Proceedings of OMAE.02 21st International Conference 
on Offshore Mechanics and Artie Engineering June 23- 
28, 2002,Oslo, Norway, number OMAE2002-28146. 

Vandemark, D., B. Chapron, J. Sun, G. H. Crescenti, 
and H. C. Graber, 2004: Ocean wave slope observations 
using radar backscatter and laser altimeters. J. Phys. 
Oceanogr., 34, 2825-2842. 

Vandemark, D., P. D. Mourad, S. A. Bailey, T. L. Craw- 
ford, C. A. Vogel, J. Sun, and B. Chapron, 2001: Mea- 
sured changes in ocean surface roughness due to at- 
mospheric boundary layer rolls. J. Geophys. Res., 106, 
4639-4654. 

Violante-Carvalho, N., F. J. Ocampo- Torres, and I. S. 
Robinson, 2004: Buoy observations of the influence of 
swell on wind waves in the open ocean. Appl. Ocean Res., 
26, 49-60. 

Wang, D. W., D. A. Mitchell, W. J. Teague, E. Jarosz, and 
M. S. Hulbert, 2005: Extreme waves under hurricane 
Ivan. Science, 309, 896. 



24 



WISE Group, 2007: Wave modelling - the state 
of the art. Progress in Oceanography, 75, 
doi:10.1016/j.pocean.2007.05.005, 603-674. 

Wu, C. H. and H. M. Ncpf, 2002: Breaking criteria and en- 
ergy losses for three-dimensional wave breaking. J. Geo- 
phys. Res., 107, doi:10.1029/2001JC001077, 3177. 

Young, I. R., 2005: Directional spectra of hur- 
ricane wind waves. J. Geophys. Res., Ill, 
doi:10.1029/2006JC003540, C08020. 

Young, I. R. and A. V. Babanin, 2006: Spectral distribu- 
tion of energy dissipation of wind-generated waves due to 
dominant wave breaking. J. Phys. Oceanogr., 36, 376- 
394. 

Young, I. R. and G. P. van Vlcddcr, 1993: A review of 
the central role of nonlinear interactions in wind-wave 
evolution. Phil. Trans. Roy. Soc. London A, 342, 505- 
524. 



25 



Parameter sec cq. variable in code WAM-Cycle4 BAJ TEST405 TEST437 TEST441 





(20) 


ALPHAO 
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0.0095 


idem 


idem 


idem 


Anax 


(21) 


BETAMAX 


1.2 


1.2 


1.55 


1.52 


idem 




(21) 


ZALP 


0.011 


0.011 


0.006 


idem 


idem 


\) , IY1 £L!JC 


(23) 
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N.A. 


N.A. 


0.002 


idem 


idem 
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idem 


idem 


Sl 


(10) 


SWELLF 
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0.8 


idem 


idem 
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(10) 
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-0.018 


idem 


idem 


S3 
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SWELLF3 
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0.015 


idem 


idem 


0.5Re c i? s 
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1 x 10 5 


idem 


idem 


idem 


Cdsv 


(8) 


SWELLF5 
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1.2 


idem 


idem 


r z o 


(11) 


ZORAT 
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0.04 


idem 


idem 




(22) 


ZWND 


10 


10 


10 


10 


10 


P 


(20) 


SINTHP 
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Table 3. Wind-wave interaction parameters as implemented in version 3.14-SHOM of the WAVEWATCH III code, and 
values used in the tests presented here. In WWATCH, all parameters are accessible via the SIN3 namclist. All of these 
parameters are included in version 3.14 of WWATCH. s is a switch that, if nonzero, activates the calculation of S out - 



Parameter 


see eq. 


variable in code 


WAM4 


BAJ 


TEST405 


TEST437 


TEST441 


TEST443 


Cds 


(2) 


SDSC1 


-4.5 


-2.1 


0.0 


idem 


idem 


idem 


r 


(2) 


WNMEANP 


-0.5 


0.5 


0.5 


idem 


idem 


idem 


/fm 


(7) 


FXFM3 


2.5 


2.5 


2.5 


9.9 


idem 


idem 


Si 


(2) 


SDSDELT 


0.5 


0.4 
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idem 


idem 


idem 


s 2 


(2) 


SDSDELTA2 


0.5 


0.6 
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idem 


idem 


idem 


°ds 


(13) 


SDSC2 


0.0 
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-2.2 x 10~ 5 


idem 


idem 


idem 


Of 


SDSLF 


1.0 


1.0 
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idem 


idem 


idem 


Cm 




SDSHF 


1.0 


1.0 


0.0 


idem 


idem 


idem 


A e 


(12) 


SDSDTH 


0.0 
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80 


idem 


idem 


idem 


B r 


(12) 


SDSBR 
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1.2 x 10~ 3 


9 x 10~ 4 


idem 


idem 


' cu 


(18) 


SDSBRF1 
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0.5 


idem 


idem 


2 * C cu 


(18) 


SDSC3 


0.0 
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-2.0 


-0.8 


idem 


SB 
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SDSCOS 


0.0 
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2.0 


idem 


idem 


s d 


(13) 


SDSDC6 
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1.0 


0.3 


0.3 


0.0 


M 




SDSBM0 
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1.0 


idem 


idem 


idem 


Mi 




SDSBM1 


0.0 
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0.2428 


idem 


idem 


idem 


M 2 




SDSBM2 


0.0 


0.0 


1.9995 


idem 


idem 


idem 


M 3 




SDSBM3 


0.0 


0.0 


-2.5709 


idem 


idem 


idem 


M 4 




SDSBM4 


0.0 


0.0 


1.3286 


idem 


idem 


idem 



Table 4. Dissipation parameter as implemented in version 3.14-SHOM of the WAVEWATCH III code, and values used 
in the tests presented here. In WWATCH, all parameters are accessible via the SDS3 namelist. The only parameters not 
defined in the present paper are C\f and which act like switches to activate the BAJ parameterization for the part 
of the spectrum with saturation below and abovce the spectrum, respectively. Most of these are also included in version 
3.14, except for C cu which is needed for the TEST437 and TEST441 with results described here. The TEST405 can be 
ran with version 3.14. The parameter Mo is a switch for the correction or not of B r into B' r , when Mq = 1, as is the case 
here, the correction is not applied. 
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